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Abstract 

In this paper we investigate expanding Bianchi type I models with two tilted fluids with linear 
equations of state. Individually the fluids have non-zero energy fluxes w.r.t. the symmetry surfaces, 
but these cancel each other because of the Codazzi constraint. Asymptotically toward the past the 
solutions approach Kasner states if the speeds of sound are less than that of light. If one of the 
fluids has a speed of sound that is less or equal to 1/3 of the speed of light (radiation) then the 
models isotropize toward the future, but if both fluids are stiffer than radiation then the final state 
is anisotropic with non-zero Hubble-normalized shear. The significance of these results is discussed 
in a broader context. 
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1 Introduction 

The construction of a relativistic model of gravity contains the following ingredients: (i) a 4-dimensional 
manifold M endowed with a Lorentzian metric, (ii) a matter source description, (iii) dynamical laws — 
Einstein's field equations and, if needed, matter equations. A general relativistic cosmological model 
requires that one in addition attempts to describe the universe at a particular scale. 

At the largest spatial scales present observational data suggest that the 'standard' ACDM model of 
cosmology provides the most simple consistent description of the universe today. This model is spatially 
homogeneous (SH) and isotropic with flat spatial geometry; the matter content consists of a dark energy 
component, modeled by a positive cosmological constant, supplemented with dark matter and atoms, 
described by pressureless fluids with n a , the unit normal to the symmetry surfaces, as the common 4- 
velocity. Density fluctuations, described by linear scalar perturbations, are seeded by an almost Gaussian, 
adiabatic, nearly scale invariant process, see e.g. Q],[2], and references therein. 

However, this description does not hold on all scales, neither spatial nor temporal. On smaller spatial 
scales matter has to be described by many components, with energy fluxes in different directions, e.g., our 
galaxy is moving w.r.t. to the CMB. In the very early universe, and perhaps also in the distant future, the 
ACDM model does not give a correct matter description, indeed, although radiation can be treated as a 
gravitational test field3 today it has been an important gravitational source in the past, and inflationary 
proponents suggest that matter could be extremely different in the very early universe. Clearly the 
matter description in the 'standard' scenario is not exact, and even with a few matter components the 
associated energy fluxes cannot all be exactly aligned, not even on the largest spatial scales. What then 
happens with the matter and its associated energy fluxes in the far future and what was the situation in 
the distant past? 

We believe that we understand how radiation and matter interacts, at least after the very early universe 
when we think our empirical experience holds. Presumably this interaction explains why the 4- velocities 
associated with radiation and matter today are fairly well aligned with each other on large spatial scales, 
or maybe this alignment was produced in the very early universe by some unknown process, perhaps 
inflation. But is it obvious that this alignment should have persisted to the extent present observations 
indicate after recombination, and is it going to persist in the far future? Do linear vector perturbations 
of SH and isotropic models suffice to determine this? In the early universe interactions presumably 
played an important role in aligning energy fluxes of different matter components, but it is unclear what 
those interactions were; is it possible to shed any light on this issue without knowing the details of these 
interactions? 

Here we are going to consider two non-interacting perfect fluids that in general move w.r.t. each other 
and a non-negative cosmological constant, which includes the ACDM matter content as a special case. 
Although this may not be a good matter description at all times, it is still a useful step since it allows a 
comparative study of the effects of various types of interactions in possible future projects, an issue we 
return to in the concluding remarks. 

Two matter source components with energy fluxes in different directions produce an anisotropic source 
which excludes the isotropic standard model and forces one to consider anisotropic geometries. There are 
several reasons that suggest that the natural anisotropic models to start with are the SH Bianchi type I 
models. One reason is their geometric simplicity since this sheds light on more general models — if type I 
turns out to be complicated, then more general models will be even worse. But more importantly is that 
they are the foundation in a hierarchy of ever more geometrically complex models. 

The SH Bianchi models (models that admit simply transitive 3-dimensional symmetry groups) form a 
crucial level in the geometric complexity hierarchy, and within this level the Bianchi type I models is the 
common ingredient since they can be obtained from all other Bianchi models by Lie algebra contractions. 
The consequences of this property are revealed when one casts Einstein's equations into a dynamical 
system where the type I models appear as part of a state space boundary that describes asymptotic 
features of all other Bianchi models, see e.g. [3] and [4], and references therein. 

Moreover, the Bianchi models themselves serve as building blocks for understanding asymptotic dy- 
namics of more general inhomogeneous models, the primary reason being the following: In the very early 
universe near a generic or isotropic singularity, or in the very late universe in an inflationary epoch, 



1 A gravitational test field is a field that does not affect the metric, i.e., it is a field for which we can neglect: (i) its 
source contribution to Einstein's equations, (ii) its influence on gravitational sources. 
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horizons form and asymptotically shrink, asymptotically prohibiting communication — a phenomenon 
naturally referred to as asymptotic silence — generically pushing inhomogeneities outside the horizons 
leading to that the asymptotic evolution can locally be described by SH models — asymptotic locality, 
see [5], [6], [7]. In the dynamical systems approach these features are formally captured by recasting the 
field equations into an infinite dimensional dynamical system that at each spatial point has a boundary — 
the silent boundary, which asymptotically attracts generic asymptotically local dynamics. The dynamics 
on the silent boundary, which determine the generic asymptotically local dynamics, is described by a 
finite dimensional dynamical system that is identical to that of the Bianchi models and hence type I is 
a common key ingredient in a very general context [5] , [6] , [7] , [9] . Furthermore, there are hints, analytical 
and numerical, that Bianchi models are important for describing future asymptotic states, even in the 
absence of inflation and asymptotic silence and locality, that sheds light on spatial structure formation. 

This is not the first study using dynamical systems techniques for studying multi-fluid models. Bianchi 
models with two fluids with both 4-velocities being orthogonal to the symmetry surfaces were studied 
in [lOj . However, it is not difficult to predict what is going to happen when the non- interacting fluids 
are aligned with each other, even in the general inhomogeneous case. Let us introduce a length scale 

1 defined by Ijl = H, where H is the Hubble scalar and where the dot refers to the time derivative 
w.r.t. the proper clock time along the common fluid congruence. For simplicity we consider fluids with 
linear equations of state such that p — wp, where p, p is the pressure and energy density, respectively, 
and where the constant w describes the speed of sound c s according to w = c 2 s when w > 0. Interesting 
examples of equations of state are: w = —1, which corresponds to a positive cosmological constant; dust, 
w = 0; radiation, w = ^; and a stiff fluid, w = 1, for which the speed of sound is equal to that of light. 
Local energy-momentum conservation then yields that p oc l~ 3 ( 1 + w ) ; see e g pj[| Thus if ^ — > the 
energy densities of fluids with smaller w become asymptotically negligible compared to those with larger 
w, while if £ — > oo the opposite holds. 

However, the situation is much more complicated when the non-interacting fluids are not aligned. 
In [12] Bianchi type V models with two fluids and a positive cosmological constant were studied; one 
of the fluids had a flow orthogonal to the symmetry surfaces while the other had a 'tilted' flow, i.e., its 
4- velocity was not aligned with the normal to the SH surfaces. In [j"3] . where brane- world cosmology 
was invoked to motivate the study of multiple fluids, Bianchi type VIo with two non-interacting tilted 
fluids was investigated. Perhaps because the focus in these papers was on the quite interesting late time 
behavior, there was no mentioning about the possibility of having two (or more) tilted fluids in Bianchi 
type I. 

The present paper is organized as follows: In the next section we derive a reduced dynamical system 
that describes Bianchi type I with two non-interacting fluids with linear equation of state and a cosmo- 
logical constant. In the subsequent section we describe the associated state space and the influence of a 
positive cosmological constant; in addition we list the invariant subsets and fix points that are essential 
for understanding the present models. In section [4] we give several monotonic functions that are useful 
for determining the asymptotic dynamics; we also briefly discuss some reasons why they exist since this 
allows one to produce monotonic functions for other models as well. Section [5] takes the results in the 
previous sections as the starting point for a dynamical systems analysis which yield our main results 
about asymptotic dynamics toward the past and toward the future, in the absence of a cosmological con- 
stant. We conclude with a discussion in section [6] about the significance of our results in a more general 
context. Appendix lAl contains detailed information about the fix points and their stability properties. 
Throughout we use units such that c = 1 = 8irG. 

2 Derivation of the dynamical system 

In the orthonormal frame approach one uses a tetrad of four orthogonal unit basis vector fields { e a } and 
the associated dual one-forms { u) a } (a = 0, 1, 2, 3), which, when expressed in a local coordinate basis, take 
the form e Q = e^d/dx^ = e a M <9 M , w a = e a M dir^ (p = 0, 1, 2, 3), where the tetrad components e Q AI (x I/ ) 
and their inverse components e a tI (x") satisfy the duality relations e Q M e a „ = <5 M „ 4=> e a M e b M = S b a , and 
where the orthogonality conditions are given by g ab = e a ■ e b = g^e^e^ = r] ab ; = r) ab e% e b l/ , 
Vab = diag ( — 1, 1, 1, 1 ). The commutator functions c a {, c (a; M ), defined as [e a ,e/,] = c c ab e c , are typically 
'elevated' to dependent variables satisfying the Jacobi identities, e[ a c d bc ] — c d e \ a c e bc ] = 0. 

Let us now consider SH Bianchi models, i.e., models with a foliation of SH hypersurfaccs invariant 
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under a simply transitive group action G3, and let us also introduce an orthonormal basis of vector 
fields {e a } that is invariant under the group action such that the timelines are orthogonal to the SH 
hypersurfaces with eo = n = d/dt, where t is the proper time along the geodesic timelines (the geodesic 
property follows from the symmetries). This yields the line-element: ds 2 = —dt 2 + 5 a p u a ® a;' 3 (a, (3 = 
1,2,3), where u) a (with components e a i) are the one-forms dual to the triad e Q (with components e a l ), 
tangential to the symmetry surfaces, i.e., e^e^i — S a @ (i = 1, 2, 3). 
A 3+1 split of the commutator equations w.r.t. eo = n yields: 

[e , e a ] = -[HdJ + a J + ej 1 W)e p , (la) 
[e Q , ep] = c 1 'a/3 e 1 = 2a [Q 5^~< + e Q/3<5 n Sl . (lb) 

Here H is the Hubble variable, which is related to the expansion 8 of n according to H = -|(9; a a p 
is the shear associated with n; f2 a is the Fermi rotation which describes how the spatial triad rotates 
with respect to a gyroscopically fixed so-called Fermi frame0 n al3 and a a describe the Lie algebra of the 
3-dimensional simply transitive Lie group and determine the spatial three-curvature, see e.g. (3J- 

Due to the symmetries e a % can be written as e a % = e a ^ ] (t)ep l , where e a l are functions of the spatial 
coordinates x l alone such that [e OJ e^] = c 7 Q/ 3e 7 = 2d^ a Sp^ + eapg n Sl , where e a = e^d/dx 1 , and where 
c 7 Q( 3, parameterized by a a , n a @ , are the structure constants of the symmetry group. The symmetries lead 
to that the equations for the variables e a ^(t) (de a ^(t)/dt = —{HS a ~ l + a a ^ + e^s Q S ]e 7 /3 , as follows 
from ([Taj) ) decouple from the remaining field equations, and because of this they are not usually considered 
when discussing Bianchi models from an orthonormal frame perspective. 

A 3+1 split of the total stress-energy tensor T a b w.r.t. n a yields: 

Tab = P n a n b + 2g (a n b ) +ph ab + ir ab , (2a) 
p = n a n b T ab , q a = -K b n c T bci p=\h ab T abl ir ab = h {a c h b) d T cd , (2b) 

where h ab = n a n b + g ab ; p,p is the total energy density and total effective pressure, respectively, measured 
in the rest space of n a ; (..} stands for the trace-free part of a symmetric spatial tensor, i.e. A^ a ^ = 

A a f3 — \8 a pA 1 1 . In general T ab consists of several components T^ b \ such that T ab — J2i^a!b ■ ^ tne 
components are non-interacting, then V a T^ = 0. A cosmological constant A can be formally regarded 



as a component of T ab such that p\ — A, pa = — A, while q" A ^ = 0, tt?R = 0. 

In the Hubble normalized approach one factors out the Hubble variable H by means of a conformal 
transformation which yields dimensionless quantities [9]. In the present SH case this amounts to the 
following: 

(£ Q(3 ,iT,A^,A a ) = jj(a a p,Q, a ,n^,a a ) , {Q,P,Q a ,Il a p) = ^(p,p,q a ,ir a p) , (3) 

where we have chosen to normalize the stress-energy quantities with 3H 2 rather than H 2 in order to 
conform with the usual definition of fl. In addition to this we choose a new dimensionless time variable 
r according to 

^ ^ w 

Since H is the only variable with dimension, its evolution equation decouples from the remaining 
equations for dimensional reasons: 

H' = -(l + q)H; q = 2E 2 + + 3P) , E 2 = \^ af3 ^ , (5) 

where a prime henceforth denotes d/dr and where q is the deceleration parameter, obtained by means 
of one of Einstein's equations — the Raychaudhuri equation; note that fl and P in the expression for q 
refers to the total Hubble-normalized stress-energy content. A 3+1 split of the remaining Einstein's field 
equations (G ab — T ab , where G ab is the Einstein tensor and T ab the total stress-energy tensor) and the 
Jacobi identities, yields the following reduced system of coupled equations for the Hubble-normalized 
variables: 



2 The sign in the definition of (1° is the same as in [8], but opposite of that in [3|. 
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= -(2 - g)S Q/ 3 + 2e 7l5 ( Q R 1 - 3 TZ( a i3) + 3£l a p > (6a) 

4, = [9 Sj - - ej\ W]A p , (6b) 

(N af3 )' = [q 5 7 (Q + 2£-/ tt + 2e 7 ( Q 5 i2*]JV^ 7 , (6c) 

0=l-S 2 + i 3 ^-fi, (6d) 

0=(3^^ + e Q ^7V <5 , 3 )S /3 7 -3Q Q , (6e) 

= Ap N^ a , (6f) 



where 3 TZ( a 0) and 3 7?. describe the trace-free and scalar parts of the Hubble-normalizcd three-curvature, 
respectively, according to: 

3 TZ (a0) =B(a/s>+2eT%JV / ,>jA r , 3 ft = -iB° a - 6^L 2 ; B a0 = 2N a7 ff - N af} ; (7) 

(|6a|) are the Hubble-normalized spatial and trace-free Einstein equations; (|6bp and (|6c| are evolution 
equations obtained from the Jacobi identities; (|6d"|) and (|6"e|) are the Hubble-normalized Gauss and Codazzi 
constraints, respectively, while the constraint |6l]) stems from the Jacobi identities. The conservation law 
V a T ah = for the total stress-energy tensor yields: 

fl' = (2q - 1) fi - 3P + 2A a Q" - £ Q/3 n Q ' 3 , (8a) 
<3« = -[2(1 - q) 5 J + + ej\ R?] Qp + (3A S a s - e af P N 7 S ) II/ . (8b) 

Let us now restrict ourselves to the Bianchi type I case with expansion (H > 0), characterized by 

A a = , N aP = . (9) 

In type I ([8b]) reduces to Q' a = —[2(1 — q) S a ^ + E a ' 5 + tcP i R 1 ] Qp, and hence the reduction of the 
Codazzi constraint (|6"e|) to Q a = is consistent since it is preserved during evolution. Thus there is no 
total energy flux in Bianchi type I, and hence the type I SH frame is an energy frame, in the nomenclature 
of Landau and Lifshitz [19j . Even so, a matter source can consist of several components that individually 
have non-zero energy fluxes, as long as they add up to zero. 

Let us now specialize to a source that consists of a non-negative cosmological constant A > and 
two non-interacting perfect fluids, i.e., Tf$ — (p^ + P(i))u?j\U^ + P(i)g ab ', V a T£x = (i = 1,2), where 
P(i),P(i): is the energy density and pressure, respectively, in the rest frame of the i:th fluid, while is 
its 4-velocity. We assume that p^) > 0, and for simplicity also a linear equations of state, pu\ = wu\p(j\, 
where wu\ = const. The most interesting equations of state are dust, w = Q, and radiation, w = but 
it is useful to not restrict oneself to these values in order to study structural stability, however, we do 
restrict ourselves to 

< w {2) < w {1) < 1 ; (10) 

since = 1, u)m = w (2) are associated with bifurcations that needs special treatment, to be dealt with 
elsewhere H 

Making a 3+1 split with respect to n = eg yields 

^) = r (i )(n° + «£)); n a vfo=0, T (l) = 1/ - (v 2 (i) = d af} v^v^) , (11) 

which gives 

of l) = (i + W(l) )(G ( ; ) )- 1 ^)%)' % = «'(i)nw + i(i-3tfl«))oLH)' n S = Q(H>' ( 12 ) 

where — 1 ± v 2 ^; il^ — p^/(3H 2 ). The cosmological constant contributes Q\ = A/(3H 2 ) — 
—Pa to the total f2 and P, while = = n^ 3 . Due to its definition and equation ([5]), Q\ satisfies 

n' A = 2(1 + q)n A . (13) 

3 We could have extended the range of the equation of state to include —1/3 < < 0, but the well-posedness of the 
Einstein equations for this range has been questioned, see 1201 . 
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The Codazzi constraint, Q a — Qa + Qa^ — 0, taken in combination with (TT2|) forces the two fluids 
3-velocities to be anti-parallel. Kinematically the situation is similar to that of Bianchi type I with 
a general magnetic field studied in 15j, and it is therefore natural to exploit the same mathematical 
structures in the present problem. We therefore choose the spatial triad so that one of the frame vectors 
is aligned with the fluid velocities, which we choose to be e^, in agreement with what is usually done 
in physics, i.e. u?k = (0,0,1^)). Demanding that these conditions on v%> hold for all times lead to the 
following conditional 

Ri = — £23 , R2 = £31 ■ (14) 

This leaves R3 undetermined, however, we still have the freedom of arbitrary rotations in the 1-2-plane, 
which we use to set 

R 3 = . (15) 
Following [15] . we introduce the variables E + , E^, E^, Ec according to 

E + = i(E 11 +E 22 ), E 31 +iE 23 = Me'*, E_ + -±=E 12 = (E B + i E c )e 2i * , (16) 

where E_ = (En — E2 2 )/(2-\/3), which leads to 

E 2 = Y? + + Y? A + J? B + E 2 . . (17) 

The above decomposition of S a( g has the advantage that the equation for 4>, dip /air — — Ec, decouples 
from the other equations, leaving the following reduced constrained dynamical system of coupled equations 
for the Hubble-normalized shear variables E_|_, S^, Es, Ec, the fluid 3-velocities W(i),i | (2) J the Hubblc- 
normalized energy densities f2( 1 ),f2( 2 )j and the Hubble-normalized cosmological constant f2 A : 

Evolution equations: 





= ~(2 


- ?)E+ + 3E 2 ! - Q(l)«(l) - Q(2)V(2) , 


(18a) 


*A 


= -(2 


— q+ 3E + + v / 3E B )E J 4 , 


(18b) 




= -(2 


-g)E B + V3Ei-2\/3E c , 


(18c) 




= -(2 


-q- 2V3E S )E C , 


(18d) 




= (G« 


, )- 1 (l-«( 4 ))(3^(. i ) -l + 2E+) W(i) , 


(18e) 




= (2?- 


- 1 - 3u;(i))n (i) + (3w(i) - 1 + 2E + )Q (l) t; (l) , 


(18f) 




= 2(1- 


r-g)0 A . 


(18g) 



Constraint equations: 

= 1 - e 2 — n^j - si( 2 ) - , 

= Q(i) + Q(2) , 

where 

q = 2E 2 + ±(Q ro + 3P m ) - fi A = 2 - §(n m - P m ) - 3ft A ; n m = + fi (2) , P m = P (1) + P (2) . (20) 

Equations (|18eP and (|18f[) were obtained by using that ((SJ takes the same form for non-interacting 
individual matter components, where, however, the total matter content enters into q, together with the 
type I conditions and the relations (|12j) for the individual perfect fluids. The assumption of non-negative 
energy densities and a non-negative cosmological constant, fi(j),f2 A > 0, together with (f20|) and (j!9a|) . 
yields that — 1 < q < 2, and hence that 2 — q > 0, where q = — 1 only when Qa = 1, f! m = 0, E 2 = 0. It 
follows that r £ (—00, 00) and H — > 00 when r — > —00 (if Qa ^= 1 initially). 

4 In the case of a magnetic field, aligned along e3, one obtains Ri = £23, R2 = — £31, i-e., the signs are opposite of 
those of the two tilted fluid case! This dynamical result in turn leads to sign differences in the S-equations. Note that the 
kinematic results in |15| still hold, hence e.g. fix points correspond to transitively self-similar models. 



(19a) 
(19b) 
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The auxiliary equation, 

P[i) = -(1 + w W )(G ( ; ) )- 1 [3 + vf l} - 2£+« 2 ) ]p (l) , (21) 

implies that p^ is a monotonically decreasing function such that p^ — > oo (pu) — * 0) when r — ► — oo 
(t — ► oo); hence the models begin with an initial curvature singularity, where A becomes negligible when 
compared to p^ when r — * — oo, and then expand forever to a state where the ordinary matter is infinitely 
diluted, leading to that fi ra becomes negligible compared to Q\- 

3 State space properties 

3.1 The state space 

The reduced state space consists of S = {E+, S^, Sc, Vm, t>(2), fim, ^12)1 ^a}j subject to the two 
constraints (|19p . i.e., the state space is seven-dimensional. From the definitions and the constraints (fT9)) it 
follows that the state space is bounded. Our primary concern in this paper is the 'interior' state space for 
the case of two tilted fluids for which f2(i)fi(2) > 0, < |u(i)W(2)| < 1, however, the solutions belonging to 
the interior state space often asymptotically approach its boundary. To understand the interior dynamics 
we therefore consider the closure of the interior state space, S, thus obtaining a compact state space, 
which is possible because of the regularity of the evolution equations. Hence £ 2 < 1,0 < v?^ < 1; 
< Hm < 1,0 < f2/2) < ljO < ft a < 1, in such a way so that the constraints (p~9|) are satisfied; note 
that the Codazzi constraint (|19b|) leads to that t>(i)U(2) < when Qmfip) > 0, a condition on that 
we extend to the boundary. The dynamical system (| 1 8|1 . (|19p is invariant under the following discrete 
symmetries: 

-* -^a , S c -> -Ec ; («(i),«(2)) -(u(i)j "(2)) • (22) 
We therefore assume without loss of generality that G [0, 1], Sc € [0, 1], vm G [0, 1], and iv 2 ) G [—1,0]; 
the solutions in the other sectors of the state space are easily obtained by means of the discrete symmetries. 

3.2 The influence of a cosmological constant 

Equation ( |18g[ ) implies that Q\ is monotonically increasing from zero to one. Hence 

fiA->l, S 2 ^0, fi w ->0 when r -> 00 , (23) 

as follows from combining Qa — > 1 with the Gauss constraint (|19a[) . i.e., the solutions approach a de 
Sitter state when t — * 00. This result is a special case of the proof by Wald [ST], which holds for Bianchi 
types I- VIII. In the present case the fluids behave as test fields on a de Sitter background at late times, 
obeying the equations: vL% = (G^) _1 (l — v ) (3w^ — 1) v^y It follows that = const if = |; 
is monotonically increasing (decreasing) from to 1 (1 to 0) if > ^ ( w (i) < f )i u (2) is monotonically 
decreasing (increasing) from to —1 (—1 to 0) if w^) > § ( w (2) < f)- Thus if one of the fluids has a 
soft equation of state, W(2) < 3, and the other has a sufficiently stiff equation of state, W(ij > ^, then the 
fluids will obtain a relative velocity w.r.t each other (in general when — | and always if > |); 
this is an invariant statement, and it is not possible to eliminate this effect with any choice of reference 
congruence — if one has two fluids, one with a sufficiently soft and one with a sufficiently stiff equation 
of state, then it follows that the fluids will asymptotically form anisotropies on a de Sitter background 
irrespectively of the choice of reference congruence. We note that this result is compatible with the 
analysis of Bianchi type V in [12] , and that it reflects a bifurcation that takes place at w = | for a fluid 
in any, homogeneous or inhomogeneous, forever expanding model with a cosmological constant, see |22j . 

At early times A has a negligible effect compared to normal matter and hence it suffices to study the 
Ha = subset (it follows from ( |18g[ ) and the application of the monotonicity principle, see e.g. [3], [T6] and 
references therein^ that the a-limit for all orbits (solutions) must reside on this subset (assuming that 

5 The monotonicity principle gives information about the global asymptotic behavior of the dynamical system. If M : 
X — » R is a C 1 function which is strictly decreasing along orbits (solutions) in X, then 

f(a:) C{(e X\X I lim Af(C) ^ sup A'/} , a(x) C {§ £ X\X | lim M(Q ^ inf M} 
<-*€ x C-»? x 

for all x £ X, where [ce(x)] is the ijj-limit [o-limit] set of a point x £ X, defined as the set of all accumulation points 

of the future [past] orbit of x; and analogously for strictly increasing monotonic functions. 
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17a 1 initially); cf. also the discussion after equation (f2Tj) ). Since A therefore has no effect on the past 
asymptotic dynamics and since it is of interest to also study late time behavior when one does not have 
a cosmological constant, we will from now on assume A = 0. The state space we henceforth therefore 
consider is given by 

S = {E + ,Ea,Ss,Sc,U(i),U(2),^(i),^(2)} , (24) 

subject to the constraints (fT9|) . i.e., the state space when one does not have a cosmological constant is 
six-dimensional; since the discrete symmetries (f22|) still hold we continue to assume that E^ S [0, 1], £c G 
[0, 1], vm € [0, 1], and «( 2 ) G [—1, 0]. When A = the deceleration parameter q is given by 

<Z = 2E 2 + i(ft m + 3P m ) -2- §(Q m -P m ) 5<9<2. (25) 

3.3 Invariant subsets 

The dynamical system (TT8l) . (fT9|l . with JIa = 0, admits a number of invariant subsets, conveniently 
divided into three classes: (i) 'geometric subsets', i.e., sets associated with conditions on the shear and 
hence the metric since the type I models are intrinsically flat; (ii) invariant sets on the boundary of the 
physical state space for two tilted fluids that do not belong to (i); (iii) subsets that can be obtained 
by intersections of the subsets belonging to (i) and (ii). We will introduce a notation where the kernel 
suggests the type of subset and where a subscript, when existent, suggests the values of and V( 2 ). 

Geometric subsets 

• TW: The 'twisting' subset, characterized by Ec = 0, E^ ^ 0, which leads to that the decoupled 
<fi- variable satisfies (f> = const and hence E12 oc En — £22- 

• 1ZT): The constantly rotated diagonal subset, given by E^ = 0, £c 7^ (R a — 0). This subset is 
the diagonal subset, discussed next, rotated with a constant angle around e%. 

• T>: The diagonal subset, defined by T,a = Ec = 0; E# = E_, and hence R a = 0. 

• C1ZS: The locally rotationally symmetric subset. This plane symmetric subset of the diagonal 
subset is characterized by the additional condition E# = E_ = 0. This is the simplest subset 
compatible with two tilted fluids. 

• !FCO, TCTq V{2) and TCT V{1) q: The demand that E 2 = 0, and hence fl m = 1, holds for all times 
enforces either that Dm = «(2) = 0, which defines the orthogonal Friedmann-Lemaitre subset TCO, 
or = $1(2) = 0, = f (U(2) = = 0, £1(2) = f), which gives the J 7 CT Qv{2) (TCT V{1)0 ) 
Friedmann-Lemaitre subset with one orthogonal fluid and a test vector field v (2) («(i)); these subsets 
belong to the boundary of the two tilted fluid case and thus there exists no Friedmann-Lemaitre 
subset with two tilted fluids. 



Boundary subsets 

• O: The orthogonal subset for which Dm = V( 2 ) = 0. In general this subset is expressed in a 
non-Fermi frame for which T,a, Ec 7^ 0, however, usually when dealing with this case one makes a 
rotation to a Fermi frame in which the shear and the metric are diagonal so that O belongs to T>. 

• OT V{1)0 and OT 0V{2) : The OT V{1)0 subset describes a single orthogonal fluid, 0( 2 ) > 0, «( 2 ) = 0, 
and a test vector field V(i) (il(i) = 0), and similarly for OTq V{2) . 

• £T 1V{2) and £T V{1)1 : The subset £T\ V{2) describes a fluid that consists of particles with zero rest 
mass moving with the speed of light = 1 => Q(i) = ^(1) = 3Pn), as follows from the Codazzi 
constraint (|19b[) ; similar statements hold for £T V 1. 

• JC: The vacuum subset is called the Kasner subset and is defined by fl m = 0; E 2 = 1; it describes 
the Kasner solutions, but in general in a non-Fermi propagated frame, and with as test fields. 

The lists above are far from complete; intersections of subsets are possible in many cases, which then 
form invariant subsets of lower dimension; an important example is: 
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• £Tn = £Ti v , 2 .n£T v . 1 .i: The double extreme tilt subset where both fluids propagate with the speed 
of light, v (1) = 1 = -U(2) (1) = 0(2) = 3P (1) = 3P (2 ), and hence = p (2) = 3p (1) = 3p (2) . 

There are also a number of fix points which we denote by a kernel that is related to a subset to which 
the fix point belong together with a subscript that indicates the fix point values of and «( 2 ); sometimes 
we also use a superscript. The fix points and their stability properties are given in Appendix [Aj here we 
give a brief summary: 

• There are a number of Kasner points, all satisfying m = 0, E 2 = 1, q = 2. The four Kasner 
circles: K° , K°, K°, K°, and the eight Kasner lines: KL± i)0 , KL± (i)1 , KLj„ (2) , KL± (2) , where 
the superscript denotes the sign of E^ . 

• There also are a number of Friedmann points with E 2 = 0. The four Friedmann points: FJq, FJ5, 
for which q — ^(1 + 3w^), and FqJ, F^J, for which q = i(l + 3w( 2 )). When W( 2 ) — | there is a 
line of fix points, FLj^ , that connects Fqq with Fi?, with q — i(l + 3w(i)), and similarly when 

(2) Z \ > 

= I then FL^ i)0 connects FqJ with F^J, with = i(l + 3w( 2 )). The superscript denotes the 
values of 0m and 0( 2 ). 

• When | < W( 2 ) < W(i) there exists two fix points: LRS„* i; i with E 2 = i(3w(i) — l) 2 , q — ^(1 + 
3w (1 )), and LRSi^ with E 2 = j(3w( 2 ) - I) 2 , g = |(1 + 3w (2 )). 

• On the twisting subset there exists the extremely tilted fix point TWn, for which S 2 = |, q = |; 
TW„< i, which exists when \ < < | with E 2 = -|(3w(i) — l)(15ui(i) — 7), q = i(l + 3w^); 

TWi u « i which exists when i < W( 2 ) < § with E 2 = j(3w( 2 ) — l)(15w(2) — 7), q = |(1 + 3w(2)). 

• Finally there exists the extremely tilted fix point Gn with E 2 = |, g = |; the line of fix points 
GL„» i exists when w^-j = | with E 2 = |t>( 1 )/(3 + 4u( 1 )), g = §, while the line of fix points GLi„* 2) 

exists when W( 2 ) = § with E 2 = f |w( 2 )|/(3 + 4|«( 2 )|), q = § . 

4 Monotone functions and their consequences 

In the analysis of Bianchi type VIo in [13] a monotone function is defined: 

= ^T^"" where %) = - ^)*(^») , (26a) 

P(2)"(2) + P(l)"(l) 

x'-|Ki)-»(2))(l-X 2 ), -1<X<1- (26b) 

The above holds whether or not we include a cosmological constant. If wm — W( 2 ), then x is a constant 
of the motion, however, here our concern is with the case > u>(2) , and then x is a monotonic function 
that increases from —1 to 1, which leads to: 

lim x = -l lim (f3/ 2 )^(2)/l3(i)^n)) = at early times (3 {2) il {2) -> , (27) 

r — ► — oo r — > — oo 

lim x = l lim (/3( 1 )f2( 1 )//3( 2 )0(2)) = => at late times 0^Q^ — * . (28) 

Combined with the Codazzi constraint (I19b[) this leads to the following possibilities if t — > — oo: 

(i) linv^-oo (O^), 0(2)) = (0,0), i.e., the solutions a-limits reside on /C, 

(ii) lim r ^_oo(0(2), vi) — (0,0), i.e., the solutions a-limits reside on OTq V{2) , 

(iii) linv^-oo V( 2 ) = —1, liniT-^-oo Q(i) = 0( 2 ), i.e., the solutions a-limits reside on £T„ (1) i, 
or combinations/intersections thereof. If r — > oo then: 

(i) lim T _^ oo (0( 1 ), 0(2)) = (0,0), i.e., the solutions cj-limits reside on JC, 

(ii) lim T ^ oo (0( 1 ) , v 2 ) — (0,0), i.e., the solutions w-limits reside on 0T„ (1) o, 



4 MONOTONE FUNCTIONS AND THEIR CONSEQUENCES 



10 



(iii) linv^oo = 1, lira,-—^ Q(i) = 0(i), i.e., the solutions cj-limits reside on ETi V{2) , 

or combinations/intersections thereof. 

Another monotonic function is given by 

V = v 2 (1) (l - vf^-^v^l - v 2 (1) )-^- w ^ ; V = 6(w (1) - w (2) )V , (29) 

where V asymptotically increases from zero to infinity. Combining V with \ to obtain a constant of the 
motion leads to Q(i)/Q(2) = const = — 1, where the latter equality is imposed by the Codazzi constraint, 
so unfortunately we obtain nothing new. However, it follows that when r — > — oo (r — > oo) then vm) — > 
or/and i>( 2 ) ~~ * — 1 ( v (i) ~ * 1 or/and W( 2 ) — ► 0), i.e., these limits also hold in the above K, cases. 

Before giving the next monotonic functions it is useful to give the following auxiliary equations: 

Q' w = 2{q - 1 + E+)Q W , 1^=2(22-1-3^))%, where (30a) 

r w = Q^(i - ^ ) ) 1 -^^ ) 2 = (i + ^)) 2 (G«)- 2 ^)(i - ■ (30b) 

If ~Ea, £c 7^ there exist two more monotonic functions: 
Mf c = Q^T^-/X d * = Q^(l - ^^-n^-iss-ss-* ; (m|»)' = 6(5 - 9w (l) )M% , (31) 

where Mj±c asymptotically increases from zero to infinity if Wi < | while it decreases from infinity to 

zero if Wi > | ; at Wi = | is a constant of the motion, reflecting that we have bifurcations when 

Wi = |, see Appendix IA1 The above four monotonic functions, x, V = T^/T^, MjQ, M^c, can be 
combined to yield three constants of the motion, but one of these is just the Codazzi constraint, so there 
only are two independent 'non-trivial' constants of the motion; here are two possible representations of 
these constants of the motion: 

C AB = (M^) 5 - 9w < 2 > (M^) 9 ™(i>~ 5 = const , D AB = V 9w w b {M^ ] c ) w ^- w ^ = const . (32) 

In addition to these monotonic functions there also exist several monotonic functions on the various 
subsets. 

The existence of monotone functions is not coincidental, a fact that will be discussed elsewhere, but 
let us here comment on Xi which is a monotonic function for all Class A models (i.e., Bianchi models 
for which A a = 0, see section [5] and e.g. [3]). Its existence is a consequence of that x is expressible 
as a dimensionless ratio of the spatial volume density and the dimensional constants £^ in class A, 
where £^ is related to particle conservation of the i:th fluid, see e.g. [23]. Interestingly there exists one 
more constant of the motion for each fluid in Class A, however, these constants of the motion, together 
with the constants only lead to an integral related to the Codazzi constraint (|19b[) . which therefore, 
unfortunately, is of no use. Incidentally, other constants of the motion exist in class B and hence, based 
on the above insight, there should exist a monotonic function also in this case, again related to particle 
conservation, but in a more complicated way. 

4.1 The /C subset 

Before continuing it is useful to discuss the Kasner subset JC. The state space of /C is given by K = 
{£ + ,£a,S_b,£c«(i),U(2)} subjected to the Gauss constraint £ 2 = Y? + + Y? A + Y? B + Yi 2 c = 1. The 
equations for the test fields vm e [0, 1], V( 2 ) S [—1,0] decouple from those of the shear and from each 
other. The state space therefore can be written as the following Cartesian product: 

K = KP x {v {1) } x {v {2) } , KP = {£+, £b, £c} , (33) 

where KP is the projected Kasner state space, which of course is subjected to £ 2 = 1. By determining 
the a- and w-limits for solutions on KP one can then determine the asymptotic states of and v ( 2 ) 
separately, and thus the a- and w-limits for solutions on /C. Let us therefore first turn to the equations 
on KP: 

£' + = 3£ 2 4 ; £^ = -(3£ + + V3Eb)E j4 ; E^ = ViT, 2 A — 2v / 3£ 2 7 ; T! c = 2\/3£b£c • (34) 
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This system admits a circle of fix points, the projected Kasner circle: KP°, see Figure 1(a) It is described 
by XU = Sc = 0, E + = E + , E_b = E_, where the constants E± satisfy E+ + E?_ = 1. 

The subset = yields that E + = E + , E^ + Eg. = E?_, where E^ is monotonically decreasing. 
The subset E c = leads to E+ - V3E B = E+ - \/3E_ while E + and E b are monotonically increasing 
and Y? A = 1 — E+ — Y? B . Projected onto the E+ — Es-plane this yields the straight lines — single frame 
transitions, using the nomenclature of [7J, given in Figures 1(b) and 1(c) (a frame transition preserves a 
Kasner state while permuting the spatial axes). As discussed in [Jj, the general case can be regarded as 
multiple frame transitions that yield the same result as combinations of single transitions which therefore 
determine the general asymptotic solution structure on JCP, for details see [7J. From this we conclude 
that the the a-limits for all solutions with E^E,^ ^ on JCP resides on the segment KP°, yielding a 



segment on KP characterized by —1 < E^ 



< 



< E_ < 



v3 
2 



i.e., the segment consists of 



sector (213) together with the fix points Q2 and T3 on KP , see Figure 1(d) 



(213) 




(132) 

(a) Kasner sectors 




11 T 2 

(b) The S^4 = subset transitions 

E B 




li ' T 2 

(c) The T,q = subset transitions 



1 ' T 2 

(d) The past attractor segment 



Figure 1: The projected Kasner circle KP° is divided into sectors (i,j,k), defined by Ej < E 3 < Ej 
where i, j, k is a permutation of 1, 2, 3, and where Ei = E + + V3E- , £ 2 = E+ — V3E_ , E3 



-2E+, and 



the points Q a , corresponding to the non-flat plane symmetric Kasner solution, and T Q , corresponding 
to the Taub form for the Minkowski spacetime. The Kasner subset is of relevance for the past dynamics 



and hence the arrows point in the past time direction in figures 1 (b) and 1 (c) which depicts single frame 
transitions projected onto E + — E^-space. The global past attractor for the general geometric set with 
E^Es ^ on JCP consists of sector (213) together with Q 2 and T 3 on KP°. 



The a-limits for solutions on fC are determined by the a-limits on JCV which determine the asymptotic 



){iwu 



1 



limits for vuy The equation for \v^\ e [0, 1] on KP° is given by: |i>(j)|' = (G_ ) 1 (1 

2E + ) \v^\. It follows that the a-limits for all orbits on JC on the general geometric set with E^Ec 7^ 
resides on the global past attractor ^4^*}, where the subscript denotes the range of values of u>(i) and 



5 FUTURE AND PAST DYNAMICS 



12 



3]} . (35a) 

i)}U{KL+ (i)1 :E+ = -!)}, (35b) 

i(3 W(1) - 1)} U {KL+ (i)1 : t + = -i(3 W(1) - 1)}U 

l),-i]}, (35c) 
i(3 W(1) - 1)} U {KL+ (i)1 : £+ = -i(3 W(1) - 1)}U 

l),-i)}U{KL+ (2) :£+ = -§}, (35d) 
i(3 W(1) - 1)} U {KL+ (i)1 : t + = -±(3u, (1) - 1)}U 

1),-1(3™ (2) -1))}U 
{KL+ (2) : £+ = -i(3 W(2) - 1)} U {Kg, : £+ G H(3w (2) - 1), -±]} . (35e) 

As toward the past, the results in [7] implies that all orbits on KfP, on the generic geometric set 
as well as all the Kasner compatible geometric subsets, asymptotically also approach KP° toward the 
future. From this it easily follows from the decoupled vu\ equations that the w-limit for any orbit on K, 
is one of the Kasner fix points on K°., KL^„. But according to the local stability analysis in Appendix [A] 
all fix points on K, are destabilized toward the future by the matter degrees of freedom in the full state 
space, leading to that the w-limit points on JC become saddles in the full state space such that no matter 
solutions with Q(i) > 0, vJ^. < 1, v%*. < 1 initially are attracted to any part of IC when r — * 00, and thus 
the w-limits for all 'interior' matter solutions either resides on OT V(1) q or £Ti V{2) such that q < 2 when 
r — ► 00, since q = 2 only on IC. 

5 Future and past dynamics 
5.1 Future dynamics 

The following theorem is easy to prove, but is nevertheless of interest. 

Theorem 5.1. If | < W( 2 ) < w (i) < 1> md if Q(i) > 0; v fi) < 1) v %) < initially, then no models 
isotropize when t — > 00, i.e., S 2 7^ when t — > 00. 

Proof. Assume that all solutions of the above type isotropize, i.e. that the w-limit set for each solution 
resides on a Friedmann-Lemaitre subset. The equations for vrn and v/ 2 ) then yield (v(i),f( 2 )) — > (1, —1) 
when t — > 00, which is a contradiction since no Friedmann-Lemaitre subset has U(i)f( 2 ) 7^ 0. Hence none 
of the solutions described in theorem 15. II isotropize when r — > 00. □ 

The above theorem does not tell us where the solutions end up when | < W( 2 ) < ^(1) < 1- This turns 
out to depend on what geometric set they belong to, leading to a division of the models into three classes: 
(i) The TZT>, T>, CR.S subsets (ii) the TW subset, and, (iii) the general case. Unfortunately we have not 
been able to prove what the global attractors are, but our local analysis in Appendix [A] together with 
numerical simulations lead to the following conjectures: 

Conjecture 5.2. If Q(x) > 0, v?-^ < 1, u 2 2 j < 1 initially, then the ui-limit for all orbits that belong to 
the geometric subsets TZT>, T>, and CTZS is the fix point LRSi„* 2) if -| < W{2) < w {i) < 1- 

Conjecture 5.3. If Q(\) > 0j v h\ < 1) v "(2) < 1 initially, then the co-limit for all orbits that belong to 
the geometric subset TW (T,a 7^ 0) is the fix point LRSi„* 2) if ^ < W(2) < \ an d ui{2) < w (i) < 1/ fi x 
point TWi„« if | < W( 2 ) < | and W( 2 ) < < 1; the fix point TWn if | < W( 2 ) < < 1. 

Conjecture 5.4. if Q(X) > 0, < 1) u ( 2 ) < 1 initially, then the ui-limit for all orbits that belong to 
the general geometric set (Ea, ^ 0) is the fix point LRSi„* 2 if \ < w (2) < 3 an d w (2) < w (i) < 1j ^ e 



U)( 2 ) , given by 



A, 2n 
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= {K?! 
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fix point TWiu» if \ < u>2 < § and u>p) < w\ < 1; ifte line of/ix points GLi„ (2) if § = wp) < W(i) < 1/ 
i/ie /ix poini Gn if | < top) < ifm < 1. 

However, models for which Q(i) > 0, v%\ < 1, t?? 2 ) < 1 initially and with < tup) < 3 do isotropize 
(this is also true if Q(i) = 0, even if the equations of state are stiffer than radiation), as shown in the 
following lemma: 



Lemma 5.5. If Q(i) > 0, i??^ < 1, v%s < 1 initially, and if < tup) 



< 



then all models isotropize 



when t 



00, i.e. 



w/ien r 



00. 



Proof. In section[3]we showed that the future w-limit of a 'matter' orbit has to reside on either OT V{1) o or 
£Tiv i2 ) with q < 2 and O m > 0. Assume that the w-limit of an orbit resides on OT V{1) q. The equations 
for the E, f2p)-variables on this subsets are just those for a single orthogonal fluid, but in general in a 
non-Fermi propagated frame. However, in a Fermi frame the single orthogonal fluid case is easily solved 

1 and E 2 — > when r — > 00. This statement is frame invariant and therefore 



and one finds that ft 



(2) 



holds for any frame, and hence it follows that the w-limit resides on the Friedmann-Lemaitre subset 
TCT v ,,,n and that E 2 — > 0. Let us now assume that the w-limit for a matter orbit resides on £T± V . , . 

1, q 



v {1) 

Then, since vm 



2 O m - 1(0(2) - P(2)) = 2 - O ro - |(1 - 3«;(2))(fi(2) - Qp)Up))> and hence 



2q - 1 - 3«J(2) =2(1-0 



'(1) ~ 2^"(2) ~ M2)J 
+ (1 - 3lW( 2 ))(l - (2) + W(2)«(2}, 



2V- "-"(2) 

> 0, since w 



(2) 



< 



where the inequality 



is strict if O m < 1, which we now assume. Then T( 2 ) in (|30b[) is strictly monotonically increasing and 
grows without bounds, but this is impossible since T( 2 ) is finite, and hence O m — > 1, and thus E 2 — > 
when r — > 00. □ 

Theorem 5.6. IfQ(i) > 0, v\\ < 1, t> 2 2 ) < 1 initially, then the u> -limit for all orbits is the fix point F[jJ 
if < u;(2) < < |; one 0/ £/ie /w; points on the line FLj} 1 if < W(i) < w (i) = /ke point F°J 

«/0 < top) < § < W(i) < 1. 

Proof. According to Lemma [53] E 2 = asymptotically toward the future. Imposing this condition on the 
£Ti V(2) subset yields the FCT V{1) q subset with um = 1, which is a special case of the other possibility 
that the w-limit of an arbitrary orbit with Qci) 7^ initially resides on the OT V{1) q subset, and hence that 
the w-limit resides on J : CT V{1) q with up) so far undetermined (0(2) = 1). To find the desired w-limit we 
only need to find the asymptotic limit of up), which, according to (|18cj) . is determined by the signature 
of 3w(i) — 1 when E 2 — 0, immediately leading to the theorem. 

The above theorems and conjectures are summarized in the global attractor bifurcation diagrams in 
figure [2l 




A 1 



A U 



(a) The TIV, V, CKS subsets 




(b) The TW subset 



(c) The general case 



Figure 2: Future global attractor bifurcation diagrams when A = for the various geometric subsets. 



5.2 Past dynamics 

Based on the local analysis in Appendix [A] the previous analysis of /C, and a numerical analysis, we make 
the following conjecture: 
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Conjecture 5.7. The a-limit for each orbit with Q(i) > 0, v?^ < 1, u5L < 1 initially on the general 
geometric set with S^Ep ^ is one of the fix points on the global past attractor A^^ for the Kasner 
subset K. given in equation (|35|) . 

Remark. For the various geometric subsets other parts of the projected Kasner circle are the relevant 
building blocks for producing the global attractor for each subset, in a similar way as in the generic case 
(e.g., in the 1ZT> case, with ~ 0, Ec 7^ 0, < E_ < 1 is the restriction on E_, in contrast to the 
generic case where < E_ < ^). 

6 Concluding remarks 

In this paper we have shown that the type I models with two tilted fluids exhibit a rich bifurcation 
structure, hinting at the complexity one can expect from models with more realistic sources and more 
general geometries. Some of our results reflect features that hold under more general circumstances, 
while others are particular for the Bianchi type I models with two-non-interacting fluids, but in this 
latter instance the present models yield a natural reference with which to compare results from more 
general settings. 

The asymptotically silent regimes of generic spacelike singularities and of an inflationary future share 
some properties: in the inflationary case all other matter fields than the inflationary field become test 
fields and do not influence the spacetime geometry — hence matter that is not inflationary matter does 
not matter for the spacetime geometry; in the case of a generic singularity fluids with speeds of sound 
less than that of light also become test fields, in this case gravity alone creates gravity to a larger extent 
than matter, and hence 'matter does not matter' in this case either [17], [TS], [5]. However, that matter 
fields asymptotically become test fields does not mean that they do not matter observationally, on the 
contrary, today to a good approximation the CMB can be regarded as a test field although it is the prime 
observational source for cosmology! 

In the present case a cosmological constant has yielded a final de Sitter state — this is a typical feature 
in a forever expanding model, as is the bifurcation at the radiation value w = ^ . Hence if one has several 
test fields, some less stiff and some as stiff or stiffer than radiation, one obtains anisotropies on a de 
Sitter background. However, one would perhaps not expect fields that are stiffer than radiation after an 
inflationary period in the early universe or in the far future, but does the bifurcation at the radiation 
value hint at that e.g. atomic matter and/or cold dark matter and radiation develop observationally 
significant relative velocities, perhaps non- linearly? As regards generic singularities, the Bianchi type I 
Kasner singularity is transformed into a singularity of 'Mixmaster' type when one considers geometrically 
more general models that admit Bianchi type II models on the silent boundary in such a combination 
with possible frame transitions so that the whole projected Kasner circle becomes unstable toward the 
past. But it is by no means uninteresting to examine the past behavior of type I, since matter sometimes 
lead to bifurcations such that matter sometimes does matter, as illustrated by e.g., a magnetic field [T5] . 
or by a kinematic description of matter [16j where matter mattered non-generically in a very subtle way, 
illustrating that it was not quite obvious that there would not be any non-generic subtle effects in the 
present case; the lack of such effects suggest that the past dynamics in general is structurally stable under 
a change from one to several fluids as long as < < 1. 

Asymptotic scenarios where non-interacting matter components do not matter may have interesting 
consequences when one introduces more realistic interacting sources. If the interactions only contribute 
source terms that are proportional to the non-interacting parts of the source, then the interactions 
presumably also become negligible for the determination of the geometry; it is only when interactions 
contribute more to the total stress-energy than the sources themselves that the matter does not matter 
property would be broken. Hence the approximation of non-interacting fields may be asymptotically less 
restrictive than one may initially think. 

There are no (quasi-) isotropic singularities, see e.g. [M],[21] and references therein, in the present 
case when E 2 ^ initially. The reason for this is that the shear completely destabilizes such singularities 
in Bianchi type I, which therefore is extremely misleading in a (quasi-) isotropic singularity context. 

The result that models with fluids stiffer than radiation asymptotically produce anisotropies toward 
the future is mathematically interesting, and shows that the isotropization results for a single fluid are 
structurally unstable within the Bianchi type I context, although from a physical point of view one would 
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not expect such equations of states at late times. The result suggests that tilted fluids may become 
as anisotropically significant as spatial curvature at late times (in the absence of inflation) when one 
considers more general models than Bianchi type I, leading to considerable complexity, further illustrated 
by the type VIo investigation in [13] , Our results about isotropization for soft equations of state may 
be regarded as a non-linear Bianchi type I generalization of perturbations of flat FRW models with two 
fluids, a reasonable approximation before dark energy has becomes significant, and it is of interest then 
to point out that one again has radiation bifurcations. 



A Fix points and local stability analysis 



In this section we use the Gauss constraint (|19a[) to eliminate f2( 2 ) globally, however, the Codazzi con- 
straint (|19b[) cannot, unfortunately, be analytically solved globally, but we can follow ch. 7 in [3] and 
use it to locally eliminate one variable, usually ^(i), at each fix point when (|19b[) is non-singular. There 
are several features that are similar for many of the fix points. All fix points, except one, have Ec = 0; 
several fix points have XU — 0. Linearization of (Il8dj) when Ec = 0, and (|18bj) when Y*a = 0, yield the 
eigenvalues 



A 5 



■[2- q - 2V3(S B )o] , \x A =-[2-q Q + 3(£+) + V^(E B )o] 



(36) 



where qo, (£b)o, (£+)o are the fix point values of q, £b, £+, respectively. For many fix points Vu> = 
or \vu)\ = 1. In these cases linearization of (|18el) yields 



A„„ 



3w {t) - 1 + 2(£+) , X 1 = -2(3w {l) - 1 + 2(E+) )/(l - w (i) ) , 



(37) 



where the subscript refers to the variable the eigenvalue is connected with while the superscript 
denotes its absolute fix point value. Let us now turn to the various individual fix points; throughout 
kernel subscripts give an indication of the absolute fix point values for vtu and v^ 2 y 

Kasner fix points: There are four circles of Kasner points and eight lines of fix points when < W( 2 ) < 



W(iy The Kasner circles are characterized by £+ = £ + , E^ = E_, E^ = Ec = 0, 
where E± are constants that satisfy E^_ + E 2 . = 1, and the following values of v^y. 



(2) 







«(i) 


= U (2) = 


0. 


K° 

^10 




= 1 , V( 2 ) = , 




K° • 


W(l) 


= , V( 2 


= -1. 


K° 
^11 


: W(i) 


= -«(2) = 1 • 


for the four 


cases 


are: 












; 


As A ; 


As c ; 




\0 . 

A D(2) ' 


3(1- 


to(i)); 3(1- 


^10 ■ 


; 


As A ; 


As c ; 


\1 . 

v (l) ' 


\0 . 

V (2) ' 


3(1- 




K° ■ 


; 


As A ; 


As c ; 




\1 . 

V (2) ' 


3(1- 






; 


As A ; 


As c ; 


A -(l) ; 


1,1 . 

V (2) ' 


2(1 + 


E+), 



' ; (2), 



where 



As, 
A?,. 



= -(3E+ + \/3E_ 



= 3w 



(i) 



1 + 2E 



+ ! 



As c 
Ai. 



= 2\/3E_ , 

= -2(3w {l) -l + 2t + )/(l-w (t) ) 



0. 

(38a) 
(38b) 

(39a) 
(39b) 
(39c) 
(39d) 



(40a) 
(40b) 



In the Kj^Q case the Codazzi constraint (119b[) is singular and hence it cannot be locally solved; in all 
other cases (|19b[) has been used to eliminate fi(i). The zero eigenvalue corresponds to that one has a 
one-parameter set of fixed points. The eight lines of Kasner fix points are characterized by E^ = Ec = 
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0, (1) = Q (2) = 0, S 2 = 1, and 



L (2) 





£+ = 


Ki 


- 370(1)) > 




= ± \ 


/l 




o < 


U(l) < 1 > "(2) = 




KI 4>1 : 


s+ = 


l(i 


- 3w(i } ) , 


s s 


= ± \ 


A 


- Y 2 


o < 


< 1 , «(2) = 


-1, 


KL o„ (2) : 


£+ = 


l(i 


- 3W( 2 )) , 




= ± \ 


A 


- T 2 


0(1) 


= 0, -1<U( 2 ) 


< o, 


KL± (2) : 


s+ = 


l(i 


- 3w (2 )) , 




= ± \ 


A 


V 2 




= 1 , -1 < U( 2 ) 


<o, 



(41a) 
(41b) 
(41c) 
(41d) 

where the superscript denotes the sign of X^. After eliminating locally by means of the Codazzi 
constraint (|19b[) . the eigenvalues for the eight Kasner lines are: 



KI 4)0 : 


0: 


0; 


; 


^s c ; 


3(1 


-w( 2 )) ; 


KL W : 


0: 


0; 


; 


^s c ; 


3(1 


~ w (i)) ; 


KI 4 (2) : 


0: 


0; 


^a ; 




3(1 


-«>(i)) ; 


KL ^( 2 , : 


0; 


0; 


^a ; 


; 


3(1 


~ w (2)) ; 



G 



-3(w ( i) - w {2) ) , 
™(i) - W(2) 



1 - W( 2 ) 
3(l0(i) - w (2) ) 



-6 



W(l) - 10(2) 



1 



W(l) 



(42a) 
(42b) 
(42c) 
(42d) 



where again = — (3E+ + \/3Sb) , As c = 2v3Xb, where X+, take the fix point values for the 
relevant line of fix points. Here one zero eigenvalue corresponds to that one has a line of fix points while 
the second is associated with the existence of a one parameter set of solutions that are anti-parallel w.r.t. 
each other on each side of the line of fix points. 

Priedmann fix points: All four Friedmann fix points satisfy E + = = = = 0, OmOp) = 
0, O m — 1, f(i)«(2) = 0. They are distinguished by their 0(j) and values according to: 

(43a) 
(43b) 
(43c) 
(43d) 



b oo ■ 


V(l) 


= 0, 


«(2) 




n (1) 


= i, 


0( 2 ) 


= o, 


plO . 




= 0, 


"(2) 


= -1, 




= i, 


^(2) 


= 0, 


^00 • 




= 0, 


«(2) 


= 0, 




= o, 


0(2) 


= 1, 


b w ■ 




= 1, 


«(2) 






= o, 


0(2) 


= 1, 



where the superscript refers to the values of On) and 0(2). The associated eigenvalues are: 



^00 • 


Al,2,3,4 = 


-1(1 


-w(i)) ; 


310(2) 


- 1 ; 


3(w ( i) -w(2)) ; 


0(i) 


eliminated , 


plO . 
i'oi ■ 


Al,2,3,4 = 


-1(1 


- ; 


310(1) 


- 1 ; 


2(1-310(2)) . 

1 - 10(2) 


n (i) 


eliminated , 










pOl . 
*00 ■ 


Al,2,3,4 = 


-1(1 


- w(2)) ; 


3w ( i) 


- 1 ; 


-3(i0(i) -w {2) ) ; 


V(2) 


eliminated , 


pOl . 

b w ■ 


Al,2,3,4 = 


-1(1 


- w(2>) ; 


3l0( 2 ) 


- 1 ; 


2(l-3i0(i)) . 
i _ ' 


o (1) 


eliminated . 



10(1) 



(44a) 
(44b) 
(44c) 
(44d) 



Here the last entry for each line of fix points refers to the variable that has been eliminated by means of 
the Codazzi constraint l|19bj) . Two of the eigenvalues of Ai^^ refer to and As c . 11 w (2) — 3 there 
exists a line of Friedmann points, parameterized by 0( 2 ), FLj° (2) , that connects FJq and FjJ- Similarly if 
10(1) = ~ there exists a line of fix points, FL^ , that connects F°J and F°q . They are given by 



FL 
FL 



10 

0f(2) 

01 

v {1) 



0(1) = , 0(2) = const , 
0(1) = const , 0(2) = , 



O(i) = 1. 
O(i) =0, 



0(2) = : 
0(2) = 1 : 



The eigenvalues associated with the two lines are: 



FL, 
FL' 



(2) 



01 



Al,2,3,4 
Al,2,3,4 



(l-W(i)); 3i0(i)-l; 0; 
(l-w( 2 )); 3i0( 2 )-l; 0; 







(i) 







(1) 



W(2) = 1/3 , 
?0(i) = 1/3 . 

eliminated , 
eliminated . 



(45a) 
(45b) 

(46a) 
(46b) 
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We now turn to fix points for which < E 2 < 1. 

Fix points on CTZS D £T V(1) i and CTZS D £Ti V(2j : When i < u>( 2 ) < W(i) there are two additional 
fix points, LRSt,* i and LRSi,;* 2 , which enter the physical state space via FjJ and FjJ when = |, 
W(2) — ^, respectively, and move into the £7?.<S-subset with increasing values of In the stiff perfect 
fluid limit (wji) = 1, W(2) = 1) the lines merge with the coalesced Kasner lines KL^ x — KL~ ' lt 
KLit, = KL7 , respectively. The two fix points are characterized by = Es = Ec = 0, and: 

(2) (2) 

LRS w * i)i: E+ = -±(3w (1) - 1) , v {1) = ^ - | , «(2) = -l, 

3(1- W(1) )(9 W(1) + 1)(1 + W(1) ) 3(1- W(1) )(5 W(1) + 1)(3 W(1) -1) 

S2 (i) = ™ i "(2) = ™ , (47a) 

v ' 32 mm v ' 32 wm 



LRSiv* : E + = -i(3i0( 2 ) - 1) , «(i) = 1 , «( 2 ) 



3W(2) - 1 
5W( 2 ) + 1 ' 



3(l-iU( 2 ))(5«;(2)+l)(3u; ( 2)-l) 3(1 - w {2) )(9w {2) + 1)(1 + w (2) ) 

— ™ ' i2 (2) — 57; • (4'D) 

v ' 32w( 2) ' 32w (2 ) 

After eliminating locally the eigenvalues for the two LRS-points are: 

LRS„ (V : A Ea = 3(2w w - 1) ; A Ss = A Sc = -|(1 - tw (1) ) ; 6 ^ 1} ~ W(2) • 



u. (2) 



(1 -«,(!)) (l± 



LRSi„ fa) : A Ex = 3(2u> (2) - 1) ; A Sb = A Sc = -§(1 - ib (2) ) ; -6 ^^ (2) 



- |(1 ~ W(2)) (l ± yjA(w {2) j) , (48a) 
where Re A{wu\) < 1; since the expression for A(u>(j)) is rather messy we will refrain from giving it. 

Fix point on TW n £T n : 

TW n : £+ = -§, E c = , E A = E B = ^ , = 1 , « {a) = -1 , fi (1) = fi (2) = ^ . (49) 

Local elimination of On) by means of the Codazzi constraint (|19b[) yields the eigenvalues: 

>--*■■ ») ; »)- « 



Fix points on TWC\£T V(1) i and TWn^T iv (2) '- When | < u?m < | there exists one more fix point on 
TW: TW„ (1) i. This fix point comes into existence when the point LRS„» 1 bifurcate into two points at 

wm — h] it then wanders away from T> when wm increases and eventually leaves the physical state space 
through TWn when wm — §■ Yet another similar fix point exists on TW if | < W( 2 ) < §: TWi«» . 
The fix points are characterized by Ec = and 

TW„ (V : E+ = -i(3 W(1) - 1) , £ A = ^§(1-«, (1) )(2«; (1) -1), Z B = V3(2w {1) - 1) , 

(! - w (i))( 1 5w(i) - 7) 
V « = W W = -25^ + 18w {1) - 1 ' V ^ = 1 ' 

n {1) =l-±{3wm-l)(15w (1) -7)-B(wm), n (2) = B(w {1) ) , (51a) 
i(3 W(2) - 1) , Eyi = v /|(l- W( 2 ) )(2u )(2) -l) , E B = V3(2 W(2) - 1) , 

(1 - W(2))(15W( 2 ) - 7) 



(2) t 2 



= 1 ' V(2) = " (2) ^ -25t«fo + l 8w(2) - 1 ' 
«(i) =5(^(2)), J2 (2) = l-i(3«; (2) -l)(15«; (2) -7) -B(itf (2 )), (51b) 
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where 

3(l-» (l) )(7-15 !C(l) )(25w[, ) -18w (l) + l) 

^ = 32(5^-5^ + 1) • (51C) 

Local elimination of yields the following eigenvalues: 

TW„ fi) i : A Sc - -|(5 - 9w (1) ) ; 6 ^~J" (2) ; A 3 , 4 ,5,6 = -f (1 - w (1) ) (l ± ^C (1) ± £> (1) ) , (52a) 

(2) 

TW lfl&) : A Ec = -|(5 - 9tu (2) ) ; -6 ^^^ ; A 3 , 4 , 5 ,6 = -|(1 - w {2) ) (l ± ^C (2) ± £> (2) ) , (52b) 

where Cu-\ — Cu\{wu\),Du-\ = D^(w^) exhibit quite messy expressions, which we therefore refrain 
from giving, such that real parts of the associated eigenvalues always are negative. 

Fix point in the generic geometric manifold: There exists one fix point Gy for which all the 
off-diagonal components of the shear are non-zero. It thus exists on the generic 'geometric' manifold, but 
on the 'matter boundary' where both fluids are extremely tilted. It is characterized by: 

G u : E + = -|, E A = i ^, S fl = S c = ^ ! « (1) =1 )V(2) = -1,0 (1) =0 (2) = 1. (53) 
Local elimination of yields the eigenvalues: 

2(9io ( 2) - 5) ^ 2(9tu ( i) - 5) 



A 1AM = -t(l ±i \/23 ± 12V5j; (54) 

At 10(2) < = | (w(2) — | < w (i)) there exists a line of fix points, GL„ (1) i (GLi v , 2 , ), connecting 
TW„. i (TWi„. ) with G n ; GL„. x and GLi„ are given by: 



. y _ 1 y _ 1 / 34tt(l)-6 1 1 / 13« (1) -6 

+_ 5 ' A ^3VfV 3 + 4, (1) ' Efi -^' Ec -3VfV 3 + 4^, 



13 < = const < 1 , V(2) = -1 



9 + 5v ( 2 1} _ 14u (1) 

n W = 3(l + t, (1) )(3 + 4t; (1) )' " (2) = 3(1 + , (1) )(3 + 4, (1) )' (55a) 



rT . V _ 1 v - 1 / - 34 ^(2) -6 y _ 1 y _ 1 /-13t>(2)-6 



3 ' 3V3Y 3-4^) ' " ^ 3V3V 3-4« 

= 1 , — 1 < «(2) = const < 13 

"W-3(l-« (2) )(3-4t, (2) )' " (2) -3(l-, (2) )(3- 4w(2) )- < "- ,ln 
Local elimination of On) yields the eigenvalues: 



(2) ovu V a-tv(2) 

6 

9 + 5v 2 . 



GL U(1)1 : ; l \ ' X ^fi = -§ (l ± y/F (1) ±G w ) , (56a) 



W (2) 

GL ll)p) : ; ^ ; A 3 , 4 ,5,6 = -| (l ± yjF {2) ± G (2) ) , (56b) 

where -FJv) = ^(^(^(j)), G(j) = Gt{\(vt{\) exhibit quite messy expressions, which we therefore refrain from 
giving, such that real parts of the associated eigenvalues always are negative. 
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